%%html
<style>
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunsx.otf');
font-weight: bold;
}
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunsi.otf');
font-style: italic, oblique;
}
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunbxo.otf');
font-weight: bold;
font-style: italic, oblique;
}
.text_cell { font-family: "Computer Modern", "Palatino", "Palatino Linotype", serif; }
.code_cell:first-child { display: none; }
h1, h2, h3, h4 { text-align: center; }
</style>
%%capture
# Set-up and formatting
%pylab inline
from scipy.integrate import odeint
import scipy.signal as signal
# picture configuration
matplotlib.rcParams["savefig.dpi"] = 300
matplotlib.rcParams["savefig.bbox"] = "tight"
from IPython.core.interactiveshell import InteractiveShell
%config InlineBackend.figure_format = "retina" # higher-def plots
Copying over necessary functions from part I.
# Constants
g = 9.8 # m/s^2
l = 0.5 # m
# Mechanics of triple pendulum system.
# Passed into odeint function.
def triple_pendulum(r, t, g, l):
phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = r # unpack variables
# precalculate trig functions
cos_1_2 = cos(phi_1 - phi_2)
cos_1_3 = cos(phi_1 - phi_3)
cos_2_3 = cos(phi_2 - phi_3)
sin_1_2 = sin(phi_1 - phi_2)
sin_1_3 = sin(phi_1 - phi_3)
sin_2_3 = sin(phi_2 - phi_3)
# define arrays and vectors
A = array([[3, 2*cos_1_2, cos_1_3],
[2*cos_1_2, 2, cos_2_3],
[cos_1_3, cos_2_3, 1 ]])
B = array([-dphi_3**2 * sin_1_3 - 2*dphi_2**2 * sin_1_2 - 3*g*sin(phi_1)/l,
-dphi_3**2 * sin_2_3 + 2*dphi_1**2 * sin_1_2 - 2*g*sin(phi_2)/l,
dphi_2**2 * sin_2_3 + dphi_1**2 * sin_1_3 - g*sin(phi_3)/l])
ddphi = matmul(inv(A), B) # angular acceleration
return concatenate(([dphi_1, dphi_2, dphi_3], ddphi)) # tack on angular velocity
def solve_triple_pendulum(r0, tmax, l=l, dt = 0.001):
t = arange(0, tmax, dt)
solution = odeint(triple_pendulum, r0, t, args = (g, l))
phi_1 = solution[:, 0]
phi_2 = solution[:, 1]
phi_3 = solution[:, 2]
dphi_1 = solution[:, 3]
dphi_2 = solution[:, 4]
dphi_3 = solution[:, 5]
return t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3
# Plots phi_1, phi_2, and phi_3 over time.
def time_plot(solution, subtitle = ""):
t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = solution
deg_1 = phi_1 * 180 / pi
deg_2 = phi_2 * 180 / pi
deg_3 = phi_3 * 180 / pi
figsize(15, 7)
fig = figure()
plot(t, deg_1, "-g", label="First pendulum")
plot(t, deg_2, "--m", label="Second pendulum")
plot(t, deg_3, ":k", label="Third pendulum")
# titles and labels
title("Pendulum Angles over Time\n{}".format(subtitle))
xlabel("Time (s)")
ylabel("Angle (deg)")
# legend
legend()
return fig
# Plots the pendulum trajectory in x-y space. NOT scaled appropriately to length l
def trajectory_plot(solution, subtitle = ""):
t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = solution
x1 = sin(phi_1)
x2 = x1 + sin(phi_2)
x3 = x2 + sin(phi_3)
y1 = -cos(phi_1)
y2 = y1 - cos(phi_2)
y3 = y2 - cos(phi_3)
figsize(15, 7)
fig = figure()
axvline(x = 0, color = "#bbbbbb")#vertical light gray line
axhline(y = 0, color = "#bbbbbb")#horizontal light gray line
plot(x1, y1, "-g", label = "First Pendulum")
plot(x2, y2, "--m", label = "Second Pendulum")
plot(x3, y3, ":k", label = "Third Pendulum")
title("Pendulum Trajectories over Time\n{}".format(subtitle))
xlabel("x / l")
ylabel("y / l")
legend()
margins(y = 0.1)
return fig
# Plots the phase space trajectories of the masses
def phase_space_plot(solution, subtitle = ""):
t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = solution
fig = figure()
title("Phase Space Trajectories of Masses\n{}\n\n".format(subtitle))
axis("off")
plot1 = fig.add_subplot(131)
axvline(x = 0, color = "#bbbbbb")# vertical light gray line
axhline(y = 0, color = "#bbbbbb")# horizontal light gray line
plot(phi_1, dphi_1, "k")
title("First Pendulum")
xlabel("Angle")
ylabel("Angular Velocity (1/s)")
plot2 = fig.add_subplot(132)
axvline(x = 0, color = "#bbbbbb")# vertical light gray line
axhline(y = 0, color = "#bbbbbb")# horizontal light gray line
plot(phi_2, dphi_2, "k")
title("Second Pendulum")
xlabel("Angle")
plot3 = fig.add_subplot(133)
axvline(x = 0, color = "#bbbbbb")# vertical light gray line
axhline(y = 0, color = "#bbbbbb")# horizontal light gray line
plot(phi_3, dphi_3, "k")
title("Third Pendulum")
xlabel("Angle")
tight_layout()
return fig
def frequency_plot(solution, max_freq = 2):
t, phi_1, phi_2, phi_3, _, _, _ = solution
# take the Fourier transforms
fft_1 = fft.rfft(phi_1)
fft_2 = fft.rfft(phi_2)
fft_3 = fft.rfft(phi_3)
# calculate our time step and use it to scale our frequency axis
dt = t[1] - t[0]
freqs = fft.rfftfreq(len(phi_1), dt)
# fft returns complex numbers, and so we take the absolute value of each
figsize(15, 7)
fig = figure()
plot(freqs, absolute(fft_1), "-g", label = "First Pendulum")
plot(freqs, absolute(fft_2), "--m", label = "Second Pendulum")
plot(freqs, absolute(fft_3), ":k", label = "Third Pendulum")
title("Fourier Transform of Pendulum Angles")
xlabel("Frequency (1/s)")
ylabel("Angle")
xlim(0, max_freq) # we don't care about high frequencies
legend()
return fig
# Generate data for contour plot
# params: the index (or indices) of the parameter in the r0 vector which to plot over.
# For example, to simulate changing the initial angle of the 2nd pendulum from 0 to pi radians, we would call
# contour_data(1, 0, pi, ...) since phi_2 is stored in index 1 of the r0 vector
# start: the starting value of the parameter
# end: the ending value of the parameter
# xSize: the number of intermediate paramter values. Will be the x-resolution
# ySize: the number of frequencies between 0 and max_freq. Will be the y-resolution
# fft_param: which parameter to perform the FFT on and to plot. Defaults to phi_3
# max_freq: Maximum frequency to save. Defaults to 2 1/s
# file: name of the file to save the data to. All data will be saved to the data/ folder
def contour_data(params, start, end, xSize, ySize, fft_param = 3, max_freq = 2, file="contour"):
initial_values = linspace(start, end, xSize)
t_max = ySize / max_freq
data = zeros((xSize, ySize))
for step, p_0 in enumerate(initial_values):
# set up initial conditions vector
r0 = zeros(6)
r0[params] = p_0
# run simulation
solution = solve_triple_pendulum(r0, tmax=t_max, dt=0.05)
fft_data = absolute(fft.rfft(solution[fft_param]))[ : ySize] # get FFT and discard high frequencies
# normalize to this solution's values, concatenate onto our results array
data[step, :] = fft_data / max(fft_data)
save("data/{}.npy".format(file), data)
def plot_contour(data, start, end, max_freq = 2):
# plot
figsize(12,7)
fig = figure()
imshow(transpose(data), extent=[start, end, 0, max_freq], origin="lower",
interpolation="spline36", aspect=0.6*(end-start)/max_freq, cmap="plasma") # 0=white, 1=black
ylabel("Frequency of Solution (1/s)")
return fig
(continued from Part I)
This last contour plot is particularly striking, as there is a very sharp discontinuity as the system moves from the periodic to the chaotic regime. It looks like this occurs at around $\phi=1.35$, but we'd like to know the exact value. We can calculate this numerically by extracting the peak value for each initial angle, and then finding where this peak value suddnely drops from around 0.45 to 0.
data_angles = load("data/contour_angles.npy")
# find peaks
peaks = [argmax(spectrum) for spectrum in data_angles]
d_peaks = diff(peaks) # find jump
drop_index = argmin(d_peaks)
# calculate angle from drop
critical_angle = (drop_index / len(data_angles)) * (3*pi/4)
critical_angle_degrees = critical_angle * 180 / pi
print("Critical angle: {:.1f}°".format(critical_angle_degrees))
Let's explore this critical angle further. Running the simulation for $\phi_1=\phi_2=\phi_3=78.8°$, we see that the motion is periodic, but barely. Notice that the phase space plots are symmetrical.
solution = solve_triple_pendulum([critical_angle, critical_angle, critical_angle, 0, 0, 0], tmax = 50)
time_plot(solution, "Critical Angle Solution").savefig("figures/critical_angle_time_sol.png")
trajectory_plot(solution, "Critical Angle Solution").savefig("figures/critical_angle_trajetory.png")
phase_space_plot(solution, "Critical Angle Solution").savefig("figures/critical_angle_phase_space.png")
frequency_plot(solution).savefig("figures/fft_critical.png");
Adding just 0.04 radians to the inital $\phi_1$, we get a chaotic solution. Now notice that the phase space plots are not symmetrical, and the Fourier plot has no strong peaks.
solution = solve_triple_pendulum([critical_angle + 0.04, critical_angle, critical_angle, 0, 0, 0], tmax = 100)
time_plot(solution, "Post-critical Angle Solution").savefig("figures/postcritical_angle_time_sol.png")
trajectory_plot(solution, "Post-critical Angle Solution").savefig("figures/postcritical_angle_trajectory.png")
phase_space_plot(solution, "Post-critical Angle Solution").savefig("figures/postcritical_angle_phase_space.png")
frequency_plot(solution).savefig("figures/fft_postcritical.png");